Dynamics of a matter-wave bright soliton in an expulsive potential 
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The stability regimes and nonlinear dynamics of bright solitons created in a harmonic potential 
which is transversely attractive and longitudinally expulsive are presented. This choice of potential 
is motivated by the recent creation of a matter-wave bright soliton from an attractive Bose-Einstein 
condensate (L. Khaykovich et at, Science 296, 1290 (2002)). The critical branches for collapse due 
to the three-dimensional character of the gas and explosion caused by the expulsive potential are 
derived based on variational studies. Particle loss from the soliton due to sudden changes in the 
trapping potential and scattering length are quantified. It is shown that higher order solitons can 
also be created in present experiments by an abrupt change of a factor of four in the scattering 
length. It is demonstrated that quantum evaporation occurs by nonlinear tunneling of particles out 
of the soliton, leading eventually to its explosion. 

PACS numbers: 05.45.Yv, 03.75.-b, 03.75.Fi 



CO 

CN 
> 

CN 

\o 
in 
o 

CN 
O 



C 

O 

o 



I. INTRODUCTION 

Solitons are a central paradigm of nonlinear physics [ 1 1 . 
They appear in systems as diverse as the ocean [2 , 
nonlinear fiber optics |3| , and Bose-Einstein condensates 
(BEC's) 0, ||. Such localized waves propagate with- 
out spreading. They retain their form after interactions 
and are robustly stable with respect to perturbation. In 
this sense they are particle-like. In fiber optics, solitons 
have become the basis of transoceanic communication 
systems [[}). They are called bright when they repre- 
sent a local maxima in the field, and dark or gray or 
kink when they represent a minima. Matter-wave bright 
solitons arise from a balance between nonlinearity and 
dispersion in the Gross-Pitaevskii equation (GPE) which 
describes the mean field of the BEC. Such a macrosopic 
bound state of atoms is intrinsically quantum mechan- 
ical, in that its existence depends on quantum, rather 
than classical motion of the atoms. Matter-wave bright 
solitons therefore represent an intriguing combination of 
nonlinear and quantum physics. 

Investigations of BEC's with attractive interactions 
have focused on their collapse [|]. In the attractive 
GPE, which for the case of constant potential is also 
known as the focusing nonlinear Schrodinger equation, it 
is well known that solutions collapse in finite time 
This collapse can be prevented, below a certain criti- 
cal nonlinearity, by the addition of a spatially varying 
external potential. Intensive studies of collapse dynam- 
ics [fl], |T^] were initially inspired by experiments on 
the attractive BEC [|l3[ |lj] and by the possibility of ef- 
ficient implementation of Feshbach resonances [[l5| . The 
latter use an external magnetic field to tune the s-wave 
scattering length which controls the strength of the non- 
linearity governing the collapse: in particular, 85 Rb has 
been the subject of recent experiments |ll], |l^]. In 7 Li 
there exists a Feshbach resonance around 725 gauss for 
the state \F — I, nip — 1) which is amenable to the 
study of stable phenomena in attractive BEC's: the scat- 



tering length for B = is positive, and becomes slightly 
negative between 150 and 525 gauss, far from the reso- 
nance and thus far from the region in which three-body 
losses are significant. 7 Li was the first species used to 
create an attractive BEC |l||, [l4| |l^|, and it is 7 Li that 
has been used to create the first matter-wave bright soli- 
tons (l]],[20|. 

The essential experimental method used in Ref. 
to transform a stable attractive BEC into a matter-wave 
bright soliton was as follows. A BEC with positive scat- 
tering length was condensed into a harmonic trap. The 
trap was adiabatically [|l] deformed into a cigar-shaped 
geometry, i.e., lu z <C u> p , where u> z and u> p are the longitu- 
dinal and axial oscillation frequencies of the atoms in the 
trapping potential, respectively. The scattering length 
was then tuned via a Feshbach resonance to be small and 
negative. The resulting condensate, a metastable macro- 
scopic quantum bound state of about 5000 particles, was 
projected onto an expulsive harmonic potential in the lon- 
gitudinal direction. It was then observed to propagate, 
without spreading, over a distance of 1.1 mm, in contrast 
to an ideal gas, which is exploded exponentially by the 
potential. Thus it was the use of the expulsive potential 
which provided clear evidence of a self-trapped, soliton- 
like state. 

In this paper the stability and dynamics of a matter- 
wave bright soliton in an expulsive potential are studied. 
The 3D GPE which describes the mean field of the BEC 
is written as 



-|lv 2 + 5 iV|^(r,t)| 2 + im 2 [^ (x 2 +y 2 ) 
d 

+u z z 2 } } tp(r, t) = ih—i/j(r, t) , 



(1) 



where g = Ai:h 2 a/m, a is the s-wave scattering length 
tuned by the Feshbach resonance, m is the atomic mass, 
N is the number of atoms in the condensate, and tjj has 
been normalized to one. Equation ([!]), in conjunction 
with variational studies, numerical integration, reduc- 
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tion to an effective one-dimensional problem, and the 
semi-classical WKB approximation, forms the basis of 
the work herein. Section || presents the size, energy, and 
critical parameters for a stable soliton. Section [nj stud- 
ies the dynamical effects particular to an expulsive poten- 
tial: the equations of motion for the soliton parameters 
are derived; weak excitations are shown to be exponen- 
tially damped; higher order solitons created by sudden 
changes in the scattering length are shown to persist in 
the presence of the potential; and an intriguing form of 
macroscopic nonlinear tunneling, which we call quantum 
evaporation, is described. The latter eventually leads to 
explosion of the soliton. A form of nonlinear tunneling 
similar to this has already been experimentally observed 
in the context of the relaxation of condensate spin do- 
mains p^] . Finally, in Sec. |y] the conclusions and ex- 
perimental outlook are presented. 



II. STABILITY REGIMES 



The first and most important question concerning the 
attractive BEC is under what conditions it remains ro- 
bustly stable. The GPE, in general, has proven a suffi- 
cient theoretical tool for describing the vast majority of 
BEC phenomena for T <C Tbec- However, as a three- 
dimensional second-order nonlinear partial differential 
equation it is not amenable to analytic methods with- 
out approximation, of which the most commonly used 
for the repulsive, or positive scattering length case, is 
the Thomas- Fermi approximation Q , which neglects the 
kinetic energy term in Eq. (Q) . In the attractive case one 
must find another method, since the kinetic energy and 
external potential terms must dominate the mean field in 
order to prevent collapse. The variational method is the 
most immediate and useful approximation p3| . A thor- 
ough variational study of the Lagrangian dynamics of a 
purely gaussian ansatz has already been made by Perez- 
Garcia et al. j24|. However, this work did not focus on 
the case of bright soliton solutions in an expulsive poten- 
tial. To this end a hyperbolic secant rather than gaussian 
ansatz shall be utilized in the longitudinal direction. In 
limits germane to experiments this choice leads to closed 
form analytic expressions for the widths and energy of 
the condensate. 

Three separate cases are therefore considered: gen- 
eral three-dimensional axisymmctric confinement, which 
is the case of Ref. lj z — 0, which allows for full an- 
alytic solutions; and finally the quasi-lD limit, which is 
the case of Ref. ||(| and the one likely to be relevant for 
future experiments on matter-wave bright solitons. 



A. Case I: Three-Dimensional Confinement 

The Gross-Pitaevskii energy functional is defined as 

Earty] = Jd 3 f jA-|V^(f)| 2 + ^|^(f)| 4 

+ ±m[u 2 p (x 2 +y 2 )+u 2 z z 2 ]\i>(r)\ 2 } • (2) 

In the case of an expulsive potential lu z < 0. A good 
stationary variational ansatz for a bright soliton solution 
is 



ex P I — ktt~ ) sech It)) ( 3 ) 



2ir Ijl 



hi 



where l p and l z are the transverse and longitudinal widths 
which form the variational parameters. Substituting 
Eq. (||) into Eq. (|J) and rescaling all variables by the 
axial oscillator frequency u) p , one finds 
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where all quantities have been scaled to the trans- 
verse harmonic oscillator frequency u> p or length o~ p = 
^h/(muj p ) : e G p = E GP /(huj p ), 7 P = l p /a p , -y z = l z /a p , 
A = Lo z /oj p , and a = Na/a p . Taking the partial deriva- 
tives with respect to 7 P and j z and assuming 7 P 7^ 1 , one 
may simplify and obtain the equations for the extrema 
points: 
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lz 



3(1 -it) 



(5) 



(6) 



Whether these extrema are saddle points, minima, 
or maxima, is determined by the Jacobean J = 
(d^Gp){dj z e GP ) - (&y p <9 7i e G p) 2 , which, using Eq. (§) 
to eliminate 7 Z , yields 



J = l 



1 



H 8 + i« 2 7p 6 + H 4 



An 2 



l0 + ^ 2 7 p 8 + (l + ^-A 2 a 4 )7, 



2 12 



(J) 



If Eq. (0) is negative for a given solution (7 po ,7zo) to 
Eqs. (|) and (||), (7 po ,7zo) is a saddle point; if it is posi- 
tive one must further examine <9 2 p £gp or d 2 ^ eqp . If these 
second partial derivatives are positive at (7 po , Izo) the ex- 
ternum is a minimum; if they are negative the extremum 
is a maximum. As Eqs. (^) and (0) are high order poly- 
nomials, it is necessary to obtain (7 po ,7zo) numerically 
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for given parameters A and a. For attractive interactions, 
such solutions are found to constitute a local minimum 
of the ecp surface for a finite range of the parameters 
and constitute bright solitons, provided that the spatial 
extent of the condensate is smaller than the longitudi- 
nal oscillator length. A specific example relevant to th e 
experiment of Ref. |L9| shall be presented in Sec. IIIB 
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Other features of interest not shown in Fig [l] are the 
critical branch for collapse for oblate traps, uj z /uj p ^> 1, 
which approaches the vertical axis as u> z /ui p — > +oo, 
and the critical point for the isotropic case lu z /lu p = 1, 
which is found by the above methods to be —Na/a p = 
0.6501 • • •. The isotropic case has been extensively stud- 
ied, and is determined from a purely gaussian hypothe- 
sis to be —Na/tjp ~ 0.671. However, the value found 
from numerical studies |26[ ^7) of Eq. ([I]) is ~ 0.575. 
The filled circles in Fig. [jjmark critical points for col- 
lapse determined by imaginary time relaxation of Eq. o), 
and differ by between 12.5% (isotropic case) and 17.5% 
(zero longitudinal potential, following section) from the 
variational predictions. A more detailed description of 
our numerical methods is made in Sec. HI. Note that 



-Na / o„ 



even lower critical values for collapse have been found in 
experiments for which the initial state differs from the 
equilibrium one ^ [l7| : the variational studies presented 
herein make similiar predictions, as shall be illustrated 



Sec. 



FIG. 1: Stability diagram for parameter space of matter- 
wave bright soliton; the upper half plane y > represents 
the case of an attractive harmonic trapping potential, while 
y < represents that of an expulsive harmonic potential. 
Solid line: explosion branch (left) and collapse branch (right) 
obtained from three-dimensional variational studies. Dashed 
line: quasi-one-dimensional variational prediction for explo- 
sion branch. Filled circles: collapse points obtained by imag- 
inary time relaxation of the Gross-Pitaevskii equation. The 
latter differ by ~ 17.5% from variational predictions. 



The critical points for collapse caused by the mean 
field, and, in the case of u z < 0, for explosion due to 
the expulsive potential, may be found by simultaneously 
solving Eq. (|) and Eq. (Q) for J = 0. It is expedient to 
use the resultant to eliminate 7^, which yields a single 
implicit equation in A and a. The numerical solution 
to this polynomial equation is plotted in Fig. [l] (solid 
line). The outer branch signifies collapse; the second, 
inner branch, which exists only for an expulsive potential, 
represents explosion. For uj z >0 the soliton can be made 
between the vertical axis and the collapse branch; for 
lu z < it can be created in the region between collapse 
and explosion. The cusp represents the doubly critical 
point beyond which a soliton cannot be formed, and can 
be calculated analytically to be 



-Na/a p 



6V2 

5 5/8 



3 - VE ~ 0.7830 . 



(a/5-2) 1 / 4 



1 



4.3113 



(8) 
(9) 



It is therefore clear that there is a limited parameter 
range under which bright solitons can be studied in 
harmonic potentials. This range, as demonstrated in 
Refs. [Ell and pi, is experimentally realizable. 



B. Case II: Longitudinal Freedom 

In the case of u z = 0, Eqs. (|J) and (0) factor to yield 

(l-7 p 4 ) 3 (7p 6 -7 p 2 + ^ 2 ) = (10) 



and J = - ( l-^)3( 7 io_ 7 6 + l a 2 7 4 + l a2)! (n) 

respectively. j p = 1, which was already discarded in 
the initial simplification leading to Eqs. Q and (^), 
corresponds to the trivial linear solution for the non- 
interacting case a = 0, since 7 P = l p /a p . The remaining 
factor in Eq. ([To|), a polynomial of order three, may be 
solved analytically to yield a local minimum of the £qp 
surface: 
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2 cos 
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/3C 



: COS 



-\-V30 



1/2 



(12) 
(13) 



where C = N 2 a 2 /a p is the dimensionality parameter, as 
explained below. The shape of the ecp surface is a bowl 
with a chute cascading down one side, a formation often 
found in high mountains: the bowl represent a stability 
basin and the chute, collapse. A second non-trivial root 
of Eq. ( p^| ) yields the saddle point at the start of the 
chute, while a third root is extraneous. The critical point 
for collapse, which is also the x-intercept in Fig. EJ, is 
contained implicitly in the argument of the arccosines in 
Eqs. @ and ©, 



Na/a p < ^74 - 0.7598 . 



(14) 
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This expression differs by only a few percent from the 
most extreme case for an expulsive potential, Eq. (||). 
The value determined numerically by imaginary time re- 
laxation of Eq. (ffl) is -Na/<jp = 0.6268 ± 0.0035, as 
illustrated in Fig.TU 

At the critical point the soliton has size (7 po ,7zo) = 
(3 -1 / 4 , 3 -1 / 4 ), from which it is clear that it is three- 
dimensional. Far from the collapse regime, i.e., for 
C < 1, Eqs. (|f) and flU) may be Taylor expanded to 
yield 
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(15) 
(16) 



For £ <C 1, 7 zo 3> 7po an d the soliton is quasi-one- 
dimensional. £ therefore characterizes the dimensionality 
of the system. 



C. Case III: The Quasi-One-Dimensional Limit 

For ( < 1 the axial harmonic degrees of freedom 
are not populated and the s ystem is essentially one- 
dimensional: 



as shown in Sec. II B , 7 p ~ 1 and j z 3> 1. 
For strong transverse harmonic confinement it is there- 
fore a reasonable approximation to assume a wavefunc- 
tion of the form ip(f,t) — <f>(z, t) exp[— (x 2 + y 2 )/(2a 2 )] 
and integrate the GPE over the transverse degrees of 
freedom. The result is called the quasi-one-dimensional 
GPE, to first order in ( 



h 2 d 2 

- ~ — -^~2 + 9idN \<p(z, t)\ 2 + huip 
2m oz z 



ih 



V(z] 
di 4 



4>{z,t) 
z,t), (17) 



where g±D = 2auj p h is the renormalized quasi-lD cou- 
pling constant |2g[] . As in Eq. ([j]), the wavefunction has 
been normalized to one. The constant factor of hu p is 
the shift in energy due to the tranverse harmonic con- 
finement, 2 x \hijjp. 

For V(z) = 0, the analytically derived criterion used 
here for the system to be in the quasi-one-dimensional 
regime is equivalent to the intuitive criterion that the 
additional energy due to atomic interactions is much less 
than the characteristic energy scale of the transverse con- 
finement, i.e., \fi — TiLj p \ <C fiujp, where /i is the chemical 
potential. Consider the general stationary solitonic solu- 
tion to Eq. (0), 



(z,t) 



1 



=sech ( — 

lz 



exp(-ifj,t/h). (18) 



Assuming the transverse portion of the 3D wavefunc- 
tion if) remains in the linear ground state, Eq. JT§| ) 
may be substituted into Eq. (|lj]) to yield \i = 



Tiujp - N 2 g 2 D m/(8h 2 



2h 2 /(m\g m \N). Then 



\(huj p /2\ < huj i 



fi - hujp 

so that C/2 < 1. 



<C fiLUp implies 
The intuitive geo- 



metrical criterion l p -C l z obtains essentially the same 
result: l p < l z = 2h 2 /(m\g m \N), so that VC/2 < 1. 
The same argument applies for V(z) ^ provided that 
its effect on the soliton length parameter is small. In 
the case where the characteristic transverse confinement 
length approaches the s-wave scattering length a, the 
GPE, Eq. (|l7j), no longer models the system and a one- 
dimensional field theory with the appropriate effective 
coupling constant must be considered instead [ p9[ . Since 
\a\ is on the order of angstroms for a stable attractive 
BEC, this regime is not relevant to the present study. 

Equation ( jig ) is a good ansatz to obtain the stationary 
properties of a solitonic solution of Eq. ( |17j ) . Consider the 
experimentally relevant case of V(z) = ^muj 2 z 2 , where 
C is retained to order one. Then, using the same varia- 
tional methods as described in Sec. II A , one obtains the 
following expression for the extrema: 



1 = 0. 



(19) 



For to 2 > there is a single real, positive root, which rep- 
resents a minimum of the energy functional; for u> 2 < 
there are two real, positive roots, which signify a mini- 
mum and a maximum. The minimum in the latter case, 



the analytic form of which shall be needed in Sec. Ill D 



written explicitly as j z = f 1 / VC where F = F(C/|A|j de- 
creases monotonically from | to 1 with increasing £/|A|, 
and is defined as follows: 



F 



G = 




(20) 

3tt 2 / 3 [1 + VT= (64tt 2 /27)(|A|/C) 2 ] 1 /3 f \\X\J 

• > X 4/3 

K 



[1 + VI - (647rV27)(|A|/C) 2 ] 1/3 ' 



r 4/3 



One obtains a critical point for explosion, either from 
the Jacobean or directly from the analytic form of the 
roots. To prevent this from occurring, the parameters 
must be chosen such that 



|ffip|A 
m\u z \/ft, 3 




hujp — C,Tiujp/2 and l z 



4.3985. (21) 



This quasi- ID prediction for stability against explo- 
sion is illustrated in Fig. [I] (dashed line). Since l z ~ 
27i 2 /(m|(7iD|A r ), this criterion implies that the longitudi- 
nal harmonic oscillator length for an expulsive potential 
must be at least abou t two times the soliton width. As 
shall be shown in Sec. Ill Dj , quantum evaporation places 
a further constraint on the available parameter space. 

In the case where v|A|/C is small, i.e., far from the 
critical value for explosion, Eq. ( p0| ) may be Taylor ex- 
panded to obtain a simple analytic expression for the 
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soliton length: 



1 + 



4 C 2 



4 C 4 



O 



c 6 



(22) 



To fourth order in y|A|/C this expansion is accurate 
to within 1% over half the permissible parameter range 
given by Eq. @. 



Non-equilibrium Soliton Evolution in Three 
Dimensions 



In order to understand the time evolution of a soli- 
ton which is created in a non-equilibrium initial state, 
one may begin with the equations of motion for the vari- 
ational parameters of Eq. (||). In general, one adds a 
quadratic imaginary factor in x, y, and z to the exponen- 
tial in Eq . (|3| ) and then performs a Lagrangian variational 
analysis pa, |24[] . Using this method, one finds 



7p 
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7p 



7 P 



2a 



(23) 



III. NONLINEAR DYNAMICS 



The dynamics of a bright soliton in an expulsive har- 
monic potential are herein studied, to the authors' knowl- 
edge, for the first time. Several intriguing new features 
arise. Firstly, sudden changes in trap parameters, such as 
was the case in Ref. , cause particle loss and shape os- 
cillations. Changing the strength of the nonlinearity via a 
Feshbach resonance has similar effects. In both cases the 
excitations are demonstrated to damp exponentially for 
small changes. Secondly, it is shown that a large change 
in the nonlinearity can be used, with present experimen- 
tal apparatus, to create a higher order soliton. Thirdly, 
the lifetime of a soliton in an expulsive potential is not 
infinite. The effective potential experienced by an atom 
in the soliton is the external potential plus the potential 
well created by the mean field, itself a function of the 
number of atoms. Thus the wavefunction undergoes tun- 
neling with a rate which depends on N. This nonlinear 
process accelerates, eventually leading to explosion. 

The exponential damping of weak excitations, nonlin- 
ear tunneling, and eventual explosion of a soliton dif- 
fer radically from the case of a constant potential. In 
the following sections these phenomena are studied an- 
alytically and numerically. In the latter case the con- 
densate wavefunction, which solves Eqs. ([!]) and ( p^ ) 
in three and one dimensions, respectively, is expanded 
as ip(r 7 t) — £j Cj(t)$j(r). The <F,(r) are used as 
a coordinate-discretized pseudo-spectral basis, and the 
Cj(t) are propagated in time with a fourth-order vari- 
able step Runge-Kutta algorithm. These numerical tech- 
niques are standard p0| . Absorbing bounds, i.e., imag- 
inary potential wells, are added far from the soliton in 
order to simulate an unbounded system, with care being 
taken to insure there are no reflections fl3l| . 

The motion of the center of mass of the condensate 
is decoupled from the relative motion of its constituent 
atoms, as is generically the case for harmonic potentials. 
The center of mass is unstable in an expulsive potential; 
in the following only the relative motion of the atoms is 
considered. Thus the condensate wavefunction is sym- 
metrized in our simulations in order to keep its center of 
mass at rest. 
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12 



3 7 2 7 2 



(24) 



where the derivative has been taken with respect to 
rescaled time t = uj p t. These equations match those of 
Ref. |24| , up to small numerical factors due to the asym- 
metry between the hyperbolic secant and gaussian func- 
tions in the ansatz. The factor of 7r 2 /12 in Eq. ( |24 ) has 
been made explicit in order to demonstrate the difference 
between the Lagrangian analysis and simply taking the 
first derivative of Eq. (|J), which gives an identical result 
up to a slight difference in the effective mass. One may 
then perform a linear perturbation analysis around the 
minimum and derive small excitation frequencies. Since 
this has already been done elsewhere Q it is not 
repeated here. 
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FIG. 2: Contour plot of energy surface for a matter- 
wave bright soliton in an expulsive potential. The param- 



eters of Ref. |L9[ have been used, with Na/a p 



-0.6750, 



the lower of the saddle points and the closed contour to the 
left within which a stable soliton may be formed; the open 
circle marks the local minimum; and the filled circle marks 
the initial conditions of Ref. Eg]: the soliton was projected 
from an attractive harmonic potential ui Z i = 2n x 100 Hz. 



The essential limitations on changes in the soliton and 
trapping parameters to avoid collapse and/or explosion 
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can be gleaned directly from the variational studies of 
Sec. [FJ. The shape of the energy surface ecpilpilz) is 
that of a long, narrow rift running along j z . The sides of 
the rift slope rapidly upwards towards infinity. Towards 
increasing "f z there is a saddle point beyond which explo- 
sion occurs; towards decreasing j z the rift turns to point 
toward the origin and ends at a second saddle point, af- 
ter which the condensate collapses. The turning is due 
to the three-dimensional nature of the soliton near col- 
lapse. Let the lower of these two saddle points have an 
energy e max . Then the curve e GP (7 p , j z ) = e max defines a 
closed loop which signifies the lower bound for creation of 
a soliton. Any soliton created with the initial condition 
7 P = and j z = which has length and width within 
this loop, though its motion along the energy surface may 
be complicated, must, by conservation of energy, be sta- 
ble against collapse and explosion. In Fig. || is shown an 
example of the energy surface for the case of Ref. Jl9| . 
The loop defined by e max is marked with a square, while 
the open and filled circles represent the minimum and 
the initial state, respectively. 

The above variational arguments do not include 
phonon or particle emission. Thus they cannot model 
the damping of excitations, which, as shall be illustrated 
in the following section is, for an expulsive harmonic po- 
tential, exponential. 



Damping of Soliton Excitations in One 
Dimension 



Bright solitons in fiber optics are well known to shed 
their excess power by radiation. This is an example of 
the robustness which makes them so useful for long dis- 
tance communications. Thus any pulse injected with an 
area of n/2 < dzT(z, 0) < 3n/2 into a fiber op- 

tic results in a soliton with amplitude oscillations which 
decay as 1/y/i, where the governing one-dimensional fo- 
cusing nonlinear Schrodinger equation is conventionally 
written as idtT = — ^d 2 T + |T| 2 T without a fixed nor- 
malization H . A sudden change in the scattering length 
or number of particles in a matter-wave bright soliton is, 
by a rescaling of the wavefunction amplitude, time, and 
position, mathematically identical to injecting a soliton 
of the form T(z, 0) = Asech(z) into the above equation, 
where A = y/rjf/rji with 



2N\a\/a 2 = \g 1D \Nm/k 2 



(25) 



from Eq. (|l7]). The subscripts i and / signify the states 
just before and just after the abrupt change in the scat- 
tering length. This case has been solved previously with 
the inverse scattering transform |53|. Supposing A = l+b 
with \b\ < i, the resulting number of particles N p emit- 
ted from a soliton with initial number of particles Ni is 



N n = N, 



Vf 



(26) 



Supposing rjf/rji = 1 + 5, where \6\ <C 1, one finds N p = 
Ni x <5 2 /4, from which it is apparent that losses are very 
small indeed. 

An identical result can be obtained by the overlap 
integral between the initial and final soliton wavefunc- 
tions jj4|, as is now demonstrated. Let (f>i(z) be the 
initial wavefunction of the soliton with parameter rji, 
and let 4>f(z) be the wavefunction of the steady state 
soliton with parameter rjf. One can always split the 
initial wavefunction <fii into a component parallel to 
(j)f and a component 5(fi± orthogonal to (fit, so that 
4>i{z) = \/l — e<j)f{z) + Sc/)±(z,0). Since e <C 1 for rjf 
close to rji, 8<p±_ has a small norm ^fl and one may 
calculate its temporal evolution by linearizing the GPE 
around the steady state solution (fit. The spectrum of 
linear excitations, also called the Bogoliubov spectrum, 
is Ek = h 2 k 2 /(2m) + where k is the wavevector 
of the emitted wave and /if is the chemical potential of 
the depleted soliton at long times |35|, [36). As a conse- 
quence of this dispersion relation, the emitted wave 5(fi± 
will spread as an unconfincd wavepacket, with a spatial 
width which scales as t and a maximum amplitude as 
Therefore the number of particles emitted by the 
soliton at long times is simply the number of particles 
in the component 6<fi±, that is, N p = Ni(l — \{cfif\(fii)\ 2 ). 
The overlap integral is 

(tfiflfa) = (<fi(z; ri f )\<fi(z; ~ 1 - S 2 /8 (27) 

to second order in S. It follows directly that N~ ~ Ni X 
(5 2 /4, as was found above by expanding Eq. (|26|) in 5. 

In Ref. the abrupt change in parameters in Eq. (|l|) 
occurs in the longitudinal trapping frequency u> z , rather 
than the strength of the mean field. Although this case 
is not covered by the inverse scattering transform above, 
one may generalize the considerations given by the over- 
lap integral of Eq. (27). One must recalculate the linear 
excitation spectrum in the presence of an expulsive har- 
monic potential. To obtain a discrete spectrum the sys- 
tem is treated as enclosed in a box of length much larger 
than the soliton size; periodic boundary conditions are 
chosen. Neglecting quantum reflections, the eigenmodes 
may be treated as semi-classical plane waves. The effect 
of the soliton on these waves is simply a shift in phase, 
similiar to the case of a free soliton |35| ], as may be seen 
by the following. 

Consider a semi-classical wave with wavevector k 
emerging from the soliton. The expulsive potential can 
be neglected over the width of the soliton provided that 
|/i| 3> fi\oJz\, i.e., |A|/C <C 1/5, so that the soliton is 
very far from explosion, as given by Eq. (|2l]). Near 
the edge of the soliton one must verify the applicabil- 
ity of the WKB approximation with the standard crite- 
rion \d/dz(k(zy 1 )\ < 1. Given k = ^/E+ 
\d/dz(k{z)- l )\ = (\u; z \ 2 z/2)/(E+\u, 
the soliton edge is being considered, z 
emerging wave is created by a global change in the soliton 
shape, k ~ l/l z so that E ~ h 2 /(2ml 2 ). This leads to the 
same condition as above, 3> 7i|a> z |. In this case the 



2 z 2 /2, 
2 z 2 /2) 3 / 2 . Since 
; L , and since the 
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transmission coefficient is T = (2kjr\ + i) 2 /(2k/r) — i) 2 — 
exp(i#) , where 9 varies from 2n to zero as 2k/rj varies 
from zero to infinity |55[ [36|. This shift in phase does 
not change the density of states in the limit where the 
box length tends to infinity. One concludes that the ex- 
pansion of linear excitations is unaffected by the soliton. 
Thus the spatial width of 5(f>± increases as exp(|w z |i) and 



its amplitude decays as exp(— \u> 
compared to l/|u; z |. Equation 



t/2) for times long as 
gives 



(28) 



In Fig. H the cases of free and expulsive harmonic po- 
tentials discussed above as well as that of an attractive 
harmonic potential are illustrated. Numerical integration 
of the quasi-lD GPE shows that a 10% change in the 
value of r\ results in irregular oscillations for the attrac- 
tive potential, damping of the maximum density as \j\fi 
for zero external potential, and damping as exp(— nt) for 
the expulsive potential, where k ~ \wz\- I* 1 the latter 
case, Fig. ||(c), N p ~ A 7 , x 0.0031, whereas the prediction 
from above is ~ 0.0014. This enhanced loss is due to 



nonlinear tunneling, as shall be explained in Sec. [II D. 
In simulations N p is determined by evaluating the norm 
squared of the wavefunction as a function of time, since 
the absorbing bounds remove atoms outside the immedi- 
ate region of the soliton. 
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FIG. 3: Fluctuations in the peak value of the quasi-lD 
wavefunction density for a 10% change in nonlinearity from 
rji = — 1.5 [im~ to rjf — — 1.35/im -1 are shown for (a) a 
constant potential, (b) an attractive harmonic potential of 
frequency u> z — 2n x 65 Hz, and (c) an expulsive harmonic 
potential of frequency u z — 2ni x 65 Hz (solid line) . (a) damps 
as 1/yi, similar to a wavepacket, (b) recurs, and (c) damps 
exponentially as exp(— nt), where k ~ |w z |: k — 72 x 2n Hz 
and 51 x 2ir Hz for the upper and lower exponential fits to 
the decaying envelope, respectively (dashed lines). The mass 
of 7 Li has been used to obtain the time scale. 



Thus, provided that the soliton is created with width 



and length dimensions that lie within the bounds ob- 
tained from variational studies of the final parameter set, 
one may expect it to shed in a time l/|w^| its excess mass 
and come to rest at the local minimum of the energy sur- 
face. 



C. Higher Order Soliton Creation 

In the preceeding section, the matter-wave parallel to 
the injection of an optical soliton of the form T(z,t — 
0) = Asech(z) with A = 1 + b and \b\ < \ was consid- 
ered. A soliton of this form is in fact the first of a denu- 
merably infinite family of higher order soliton solutions 
to the one-dimensional focusing nonlinear Schrodinger 
equation. Each such solution can be produced from the 
initial state A = n, where n is a positive integer, with- 
out loss of particles. Moreover, provided that \b\ < |, 
A = n + b will always result in a soliton of order n, 
with a fractional loss of N p /Ni = b 2 j(n + b) 2 , as de- 
rived in Ref . . Bright solitons with zero relative phase 
have an attractive interaction. As a soliton of the form 
T(z, t = 0) = 7isech(;z) has less than n times the energy 
of a soliton of form T(z,t = 0) = sech(z), a higher or- 
der soliton can be interpreted as tightly bound state of 
n strongly overlapping solitons. It is periodic in time in 
both its phase and density. 

A matter-wave higher order soliton can be created by 
a sudden shift in the scattering length. In particular, 
the order two soliton requires that the scattering length 
be decreased by a factor of four: for 7 Li this factor can 
be produced simply for the state \F = l,mp = 1) by a 
shift in magnetic field from ~ 175 gauss to ~ 275 gauss, 
thereby shifting a from — ao to — 4ao, where ao is the 
Bohr radius. As with the soliton of order one, higher 
order solitons persist in the presence of a sufficiently weak 
expulsive potential. 

One supposes two experimental approaches: in the 
first, a soliton, already formed in an expulsive potential, 
is subjected to a sudden change in the scattering length; 
alternatively, the scattering length of an attractive con- 
densate trapped in an attractive harmonic potential is 
simultaneously shifted by a factor of four with the pro- 
jection of the soliton onto a weakly expulsive potential. 
The latter may be easier for observing oscillations of the 
soliton shape, since the center of mass in unstable unless 
the relative position of the potential and the soliton are 
continuously adjusted. In either case the initial scatter- 
ing length must be chosen such that it satisfies Eq. (pi]), 
lest the soliton explode as it returns periodically to the 
form of its initial state. Thus, to avoid both collapse and 
explosion on experimental time scales and simultaneously 
accomodate the factor of four change in 77, one should 
choose trapping potentials such that \uj z /lu p \ 1 / 2 < 0.1, as 
deduced from Fig. 

The order two soliton solution to Eq. ( |l7| ) with V(z) = 
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takes the form fl3J 



(z,t) 



exp 



128m 



3ex P [ l i6k 



+ cosh 2fii 



cosh 



V 8 



3 cos 



^ 16m J 



4 cosh 



cosh 



"Mi 

v 2 



(29) 



where 77 ^ = 2A r |aj|/cr 2 ,, as above, with a,f being the final 
scattering length. Thus the period is 16m / ^irUr) 2 ) while 
the spatial width is on the order of S/rjf. Figure [| illus- 
trates the creation of an order two soliton by abruptly 
changing the scattering length from r\i = — 0.1875 ^m -1 
to r\j — — 0.75 /mi -1 and simultaneously projecting the 
soliton from a constant to an expulsive potential of fre- 
quency oj z — 27r i x 5 Hz. This matches the analytic form 
of Eq. (P9|) . A negligible number of particles N p /Ni ~ 
0.001% are emitted due to the projection, in agreement 
with the prediction of Eq. (p8| ) , here extended to the case 
n = 2. The above simulated parameters are equivalent to 
an experiment for which uj p = 1000 x 2n Hz, uj z = 2iri x 5 
Hz, dj = — ao, a,f = — 4ao, and N — 2550 particles for 
7 Li. The shape variations of this time-dependent nonlin- 
ear structure should be observable in situ with a period 
of 19 ms. 



D. Quantum Evaporation 



l z = 2F/r] by the real positive root of Eq. ( |Tgj ) for u> 2 < 0, 
and the chemical potential by subsequently minimizing 
the eigenvalue of Eq. (|l7|). The semi-classical WKB tun- 
neling rate, scaled to the frequency of the expulsive po- 
tential, is then given by 



r 



c 



i/'i 



F|A| h\u, 



(a) 



exp(-2J) . 



(31) 
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A bright soliton in an expulsive potential experiences 
an effective potential U(z, t) — mujlz 2 /2+gi£>N\(f>(z,t)\ 2 , 
as sketched in Fig. [|. This is a time-dependent potential 
well, the depth of which is a function of the number of 
particles. The soliton therefore undergoes nonlinear tun- 
neling, which we have called quantum evaporation. To 
estimate the rate of quantum evaporation one may con- 
sider linear tunneling in the external potential U(z,t), 
where one neglects the explicit time dependence of the 
wavefunction. The semi-classical WKB approximation 
then predicts, for a one-dimensional system [B7J, 



r = 2- exp -2 



dzp(z) 



(30) 



where the first factor of 2 comes from the fact that the 
wavefunction can tunnel in two directions, a; is a charac- 
teristic oscillation frequency of the particles in the well, 
b — a is the barrier length through which the wavefunc- 
tion must tunnel, and p(z) = \j2m^J fi — U(z), where [i is 
the chemical potential. Since the effective potential well 
contains a single bound state (see Fig. ||), one may take 
heuristically \ fiuj — |?7|</>(0)| 2 — n\, which is the height of 
the chemical potential above the minimum of the well. 

Supposing that the wavefunction is modeled by the 
variational ansatz Eq. (|l8|), one may determine fi and 
the well depth 77|0(O)| 2 directly from the methods out- 
lined in Sec. II C[ 4>{z) is given by Eq. (fL8|), the length 



00 
H 
< 



-71/2 




-32 



16 



16 32 

z ([i m) 



FIG. 4: The ID numerical evolution of the (a) density 
and (b) phase of a higher order soliton created by an abrupt 
change in scattering length by a factor of four. The oscilla- 
tion period is 19 ms for 7 Li, with the initial state (solid line), 
quarter period (dotted line), and half period (dashed line) 
shown. An expulsive potential hill of frequency ui z — 2ni x 5 
Hz is present in the simulation, with r]i — 0.187 /im _1 and 
77/ = 0.750 /im _1 . Supposing a transverse oscillator frequency 
of Lu p = 2n x 1000 Hz and a change in scattering length from 
at = — ao to 0/ = — 4ao, the resulting soliton consists of 2550 
atoms and is well within the quasi-one-dimensional regime. 
Note the phase jump of n in the phase corresponding with 
the spontaneous development of nodes in the density. Such 
a nonlinear time-dependent structure could be created and 
observed in present BEC experiments. 
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where 



I = I dx 



2F 2 |A| M , F 4 |A| 2 2 



x 2 + 2Fsech (x) 



and 



C%«l c 2 

C 2C tt 2 |A|F 2 



1/2 



% z | V6|A|F 2 3|A|F 24( 



(32) 
(33) 



Here the integration variable x = z/^ and the parame- 
ters a and b are determined by U(a), U(b) = /z, as illus- 
trated in Fig. §. The function F = F(C/|A|), as defined 
by Eq. (p(j|), varies monotonically between 1 and 4/3, 
the latter value occuring at the critical point for explo- 
sion. Here the additional chemical potential tko p due to 
transverse confinement is neglected, since in the quasi-lD 
regime this is assumed to have no effect on the dynam- 
ics. Equation (|3 1|) is a function of the single parameter 
f/|A| = N 2 a 2 upn/(h\u z \). 

In the case where the soliton is far from explosion, i.e., 
C/|A| » v/267r 2 /3 3 = (C/|A|) C , as given by Eq. ©, the 
above reduces to the following simple limit: F — > 1 + , 
li/(h\u a \) - — C/(2|Ai), and / -» </(4|A|) - V21n(l + 
•\/2) + 0(|A|/C). The asymptotic tunneling rate is there- 
fore 



r 



2(1 + V2) 2V ^ 



I Ml 



exp 



I Ml 



(34) 



a function of the single parameter /i/(7i|cj z |). 

In Fig. |^(a) is shown the complete and asymptotic 
WKB variational predictions of Eq. (|3l"|) (solid line) and 




5.0 10.0 
z (|im) 



FIG. 5: The effective potential well experienced by a bright 
soliton in a one-dimensional expulsive harmonic potential is 
illustrated. The energy of the single bound state is fj,, and the 
wavefunction tunnels in either direction through a barrier of 
length b — a, as indicated by the dashed lines. Here we have 
used the parameters r\ = — 3.5/im" 1 with a steep potential 
hill of ui z = 2iri x 450 Hz for illustration purposes. 



Eq. (|34|) (dashed line), respectively, and numerical data 
obtained by direct simulation of Eq. ([l7]) (black dia- 
monds). In Fig. ^(b) the nonlinear tunneling of a soli- 
ton in an expulsive harmonic potential is shown explic- 
itly. The solid line shows the result given by numeri- 
cal evolution of the qu asi-lD G PE. In this simulation 
(C/|A|)i = 7.1151, or 2y/(C/\X\)i = 5.3348. It is apparent 
that the tunneling accelerates and the soliton explodes, 
the latter occurring after about 5% of the particles in 
the soliton have been lost. This differs by about 10% 
from the predicted critical point of Eq. (|2l|), pointing to 
a level of accuracy in the quasi- ID variational approxima- 
tion si miliar to that of the 3D variational approximation 
of Sec. II A and Fig. ||. The dashed line represents the 
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FIG. 6: Quantum evaporation, (a) The WKB prediction 
for the tunneling rate as a function of £/|A| = N 2 a 2 \uj z \m/h 
(solid line), and its large |/i|/(7i|a)z|) asymptotic expression 
(dashed line). Black diamonds: numerical data produced by 
integration of the quasi-lD GPE. (b) Solid line: explosion 
of a matter-wave bright soliton in the quasi-one-dimensional 
regime near the critical point, obtained by numerical inte- 
gration. Dashed line: WKB prediction of the same. Here 
rji = — 3.0/im _1 and uj z — 2ivi x 450 Hz. The mass of 7 Li has 
been used to obtain the time scale. 



10 



numerical integration of dN/dt = —TN, where T is ob- 
tained by the WKB approximation as shown in Fig. [|(a). 
The explosion time is similar to that of the simulation; 
the critical particle number is as given by the variational 
approximation, as assumed in the above derivation. 

To experimentally observe quantum evaporation and 
subsequent explosion one faces the difficulty that the cen- 
ter of mass of the soliton is unstable. There are two solu- 
tions: one may select parameters such that the tunneling 
rate is maximized and try to observe the soliton as it 
falls; or one may use a soliton which is large enough to 
observe in situ and construct a feedback loop to contin- 
uously shift the peak of the expulsive potential so as to 
keep the soliton centralized at its maximum. 

Note that a soliton of order n > 1 is expected to show 
the same kind of tunneling properties. 



IV. CONCLUSIONS 

The size, energy, and stability regimes of a stable 
matter- wave bright soliton in a harmonic potential which 
is transversely attractive and longitudinally expulsive 
have been presented. The critical branches for collapse 
and explosion have been delineated analytically. The 
equations of motion of the soliton parameters were calcu- 
lated, and the response of the soliton to sudden changes 
in trap parameters or the scattering length was shown to 
result in excitations which damp exponentially by emis- 
sion of particles from the soliton. It was shown that a 
higher order soliton can be created by a sudden change 
in the scattering length by a factor of four. Finally, an 
intriguing nonlinear dynamical effect was presented: a 
matter-wave bright soliton undergoes quantum evapora- 
tion and eventually explodes. 



The parameter space for creation of a locally stable 
matter-wave bright soliton is accessible to present exper- 
iments, as already demonstrated in Refs. [jl9| and [ p0[ . 
Data from the former shows that solitons are self- 
adjusting: when the scattering length of a condensate 
with ~20,000 atoms was tuned to be negative, a large 
thermal cloud with a stable condensate core of ~5,000 
atoms was observed to remain. Projection onto an ex- 
pulsive potential caused the observed thermal cloud to 
be rapidly expelled from the vicinity of the condensate, 
resulting in a nearly pure soliton. The choice of a weaker 
expulsive potential and a less negative scattering length 
than that of Ref. [|l9| could be used to create solitons 
observable in situ. Thus the experimental outlook for 
the investigation of matter-wave bright solitons is very 
positive. 

An important aspect of solitons which could be investi- 
gated in future work is their interactions: this is an essen- 
tial aspect of solitons which gives them their particle-like 
nature, in contrast to a solitary wave. Binary interac- 
tions between solitons depend on their relative phase: 
for A0 = they are attractive, while for A</> = tt they 
are repulsive |38|, [39|. Such interactions could therefore 
measure phase decoherence times. 
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